/****************************************************************
Compile: gcc -O3 -o pararobust2.c pararobust2.c
Usage: ./pararobust2 mode mode1 mode2 set pnflag N L inputfile
       ./pararobust2 0 0 0 1 0 10 13 Budding2.txt

mode          : 0 normal ; 1 cycled ; 2 triggled cycled
mode1         : 0 disallow self-activation ; 1 allow self_activation
mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II
set           : 1-5 for different parameter set
pnflag        : 1 print the solutions
N             : number of the nodes
L             : number of the steps        

This program is used to calculate the designability for a given trajectory
under the threshold rules.
****************************************************************/

#include <stdio.h>
#include <stdlib.h>

int *robustdist;
double *robustDist;

int minlink, minlinks;

struct parameter {
	int Rij, Rii, Gij, Gii, theta;
} para;

void print_help ()
{
	printf ("Usage         : ./pararobust2 mode mode1 mode2 set pnflag N L inputfile\n");
	printf ("                (./pararobust2 0 0 0 1 0 10 13 Budding2.txt)\n");
	printf ("mode          : 0 normal ; 1 cycled ; 2 triggled cycled\n");
	printf ("mode1         : 0 disallow self-activation ; 1 allow self_activation\n");
	printf ("mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II\n");
	printf ("set           : 1-5 for different parameter set\n");
	printf ("pnflag        : 1 print the solutions\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n\n");
}

void calc_links (int N, int solution[N])
{
	int i, links;

	links = 0;
	for (i = 0; i < N; i++)
	{
		if (solution[i] != 0)
			links ++;
	}
	if (links < minlink)
		minlink = links;
}

void print_solution (int N, int solution[N], int i)
{
	int j;

	printf ("Node %2d :", i);
	for (j = 0; j < N; j++)
	{
		if (solution[j] > 0)
			printf (" g");
		else if (solution[j] < 0)
			printf (" r");
		else
			printf (" n");
	}
	printf ("\n");
}

int check_sum (int N, int L, int state[L][2 * N + 1], int solution[N], int i, int k, int flag, int theta)
{
	int l, sum;

	sum = 0;

	if (flag == 0)
	{
		for (l = 0; l < N; l++)
		{
			if (l == i)
				sum += solution[l];
			else
				sum += state[k][l] * solution[l];
		}
	}
	else
	{
		for (l = 0; l < N; l++)
		{
			if (solution[l] == 0)
				continue;
			else if (i == l)
				sum += solution[l];
			else if ((solution[l] > 0 && state[k][N + i] == 1) || (solution[l] < 0 && state[k][N + i] == 0))
			{
				sum += state[k][l] * state[k][N + l] * solution[l];
			}
			else if (state[k][l] + state[k][N + l] != 0)
			{
				sum += solution[l];
			}
		}
	}

	if ((sum > theta && state[k][N + i] == 1) || (sum < theta && state[k][N + i] == 0) || (sum == theta && state[k][N + i] == state[k][i])) 
		return 1;
	else
		return 0;
}

int check_sum1 (int N, int L, int state[L][2 * N + 1], int solution[N], int i)
{
	int k, l, sum, flag[N];

	for (k = 0; k < N; k++)
		flag[k] = 0;

	for (k = 0; k < L; k++)
	{
		for (l = 0; l < N; l++)
		{
			if (l == i)
				continue;
			if (state[k][l] * solution[l] != 0)
			{
				if ((state[k][N + i] == 1 && solution[l] < 0) || (state[k][N + i] == 0 && solution[l] > 0))
				{
					flag[l] = 1;
				}
			}
		}
	}

	sum = 0;

	for (k = 0; k < N; k++)
		sum += flag[k];

	return sum;
}

	/*
		for (l = 0; l < N; l++)
		{
			else if ((solution[l] > 0 && state[k][N + i] == 1) || (solution[l] < 0 && state[k][N + i] == 0))
				sum += state[k][l] * state[k][N + l] * solution[l];
			else if (state[k][l] + state[k][N + l] != 0)
				sum += solution[l];
		}
	*/

int calc_design2 (int N, int L, int mode, int mode1, int mode2, int pnflag, int state[L][2 * N + 1], int i, int solution[N], int j)
{
	int k, l, NS, sum, flag;

	if (j == N)
	{
		for (k = 0; k < L; k++)
		{
			if (mode == 2 && k == 0)
				flag = 0;
			else if (mode2 == 0)
				flag = 0;
			else if (mode2 == 1)
			{
				if (state[k][N + i] == state[k][i])
					flag = 0;
				else
					flag = 1;
			}
			else if (mode2 == 2)
			{
				if (state[k][N + i] != state[k][2 * N])
					flag = 0;
				else
					flag = 1;
			}
			else if (mode2 == 3)
			{
				if (state[k][N + i] == state[k][i] && state[k][N + i] != state[k][2 * N])
					flag = 0;
				else
					flag = 1;
			}

			if (check_sum (N, L, state, solution, i, k, flag, para.theta) == 0)
				return 0;
		}

		sum = check_sum1 (N, L, state, solution, i);
		robustdist[i * (L + 1) + sum] ++;
		calc_links (N, solution);
		if (pnflag)
			print_solution (N, solution, i);

		return 1;
	}

	NS = 0;
	if (i != j)
	{
		solution[j] = para.Rij;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		solution[j] = 0;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		solution[j] = para.Gij;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
	}
	else
	{
		solution[j] = para.Rii;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		solution[j] = 0;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		if (mode1 == 1)
		{
			solution[j] = para.Gii;
			NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		}
	}

	return NS;
}

int calc_design (int N, int L, int mode, int mode1, int mode2, int pnflag, int pathway[L][N], int i)
{
	int j, k, state[L - 1][2 * N + 1], solution[N];

	for (j = 0; j < L - 1; j++)
	{
		for (k = 0; k < N; k++)
		{
			state[j][k] = pathway[j][k];
			state[j][N + k] = pathway[j + 1][k];
		}

		if (j != L - 2)
			state[j][2 * N] = pathway[j + 2][i];
		else
		{
			if (mode == 1)
				state[j][2 * N] = pathway[1][i];
			else
				state[j][2 * N] = pathway[j + 1][i];
		}
	}

	if (mode == 2)
	{
		if (i == 0)
			return 1;
		else
		{
			if (pathway[0][i] == pathway[1][i])
				solution[0] = 0;
			else if (pathway[1][i] == 1)
				solution[0] = para.Gij;
			else if (pathway[1][i] == 0)
				solution[0] = para.Rij;

			return calc_design2 (N, L - 1, mode, mode1, mode2, pnflag, state, i, solution, 1);
		}
	}
	else
	{
		return calc_design2 (N, L - 1, mode, mode1, mode2, pnflag, state, i, solution, 0);
	}
}

void combine_dist (int N1, int N2, int dist[N2 + 1], double Dist[N1 * N2 + 1], int i)
{
	int j, k;
	double Dist2[N1 * N2 + 1];

	if (i == 0)
	{
		for (j = 0; j <= N1 * N2; j++)
			Dist[j] = 0.0;

		for (j = 0; j <= N2; j++)
			Dist[j] = dist[j];
	}
	else
	{
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist2[j] = 0.0;
			if (j <= N2)
			{
				for (k = 0; k <= j; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
			else
			{
				for (k = 0; k <= N2; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
		}

		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist[j] = Dist2[j];
		}
	}
}

void doprint_stat (int M, double Dist[M])
{
	int i, min, max, peak, median;
	double avg, sum, sum2, maxdist; 

	min = M; max = 0;
	avg = 0.0; sum = 0.0; maxdist = 0;
	for (i = 0; i < M; i++)
	{
		if (Dist[i] > 0.1)
		{
			if (i < min)
				min = i;
			if (i > max)
				max = i;
			if (Dist[i] > maxdist)
			{
				maxdist = Dist[i];
				peak = i;
			}
			avg += i * Dist[i];
			sum += Dist[i];
		}
	}

	sum2 = 0.0;
	for (i = 0; i < M; i++)
	{
		sum2 += Dist[i];
		if (sum2 * 2 > sum)
		{
			median = i;
			break;
		}
	}

	avg /= sum;
	avg /= M;
	printf ("c %3d %3d %3d %3d %3d %.3lf ", min, max, median, peak, (int) (avg * M + 0.5), avg);
}

int main (int argc, char *argv[])
{
	FILE *fp;
	char buf[256];
	int i, j, N, L;
	int mode, mode1, mode2;
	int pnflag, parameter;

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	mode = atoi (argv[1]);
	mode1 = atoi (argv[2]);
	mode2 = atoi (argv[3]);
	parameter = atoi (argv[4]);
	pnflag = atoi (argv[5]);
	N = atoi (argv[6]);
	L = atoi (argv[7]);

	robustdist = (int *) calloc (N * L, sizeof (int));
	robustDist = (double *) calloc (N * (L - 1) + 1, sizeof (double));

	if (parameter == 1) /* strong inhibition rule */
	{
		para.Rij = -2 * N;
		para.Rii = -1;
		para.Gij = 2;
		para.Gii = 1;
		para.theta = 0;
	}
	else if (parameter == 2) /* majority rule */
	{
		para.Rij = -2;
		para.Rii = -1;
		para.Gij = 2;
		para.Gii = 1;
		para.theta = 0;
	}
	else if (parameter == 3)
	{
		para.Rij = -2;
		para.Rii = -2;
		para.Gij = 2;
		para.Gii = 2;
		para.theta = 0;
	}
	else if (parameter == 4)
	{
		para.Rij = -2;
		para.Rii = -2;
		para.Gij = 2;
		para.Gii = 2;
		para.theta = 1;
	}
	else if (parameter == 5)
	{
		para.Rij = -2;
		para.Rii = -2;
		para.Gij = 2;
		para.Gii = 2;
		para.theta = -1;
	}
	else
	{
		printf ("Unknown parameter set!\n");
		exit (-1);
	}

	int pathway[L][N], NS2[N], NS2_min;
	double NS;

	fp = fopen (argv[8], "r");
	while (fgets (buf, 256, fp))
	{
		for (i = 0; i < L; i++)
		{
			for (j = 0; j < N; j++)
			{
				fscanf (fp, "%d", &pathway[i][j]);
			}
			fgets (buf, 256, fp);
		}

		for (i = 0; i < N * L; i++)
			robustdist[i] = 0;

		NS = 1.0; minlinks = 0;
		for (i = 0; i < N; i++)
		{
			minlink = N;
			NS2[i] = calc_design (N, L, mode, mode1, mode2, pnflag, pathway, i);
			NS *= NS2[i];
			if (NS2[i] == 0)
				break;
			if (i == 0 || NS2[i] < NS2_min)
				NS2_min = NS2[i];
			minlinks += minlink;
		}

		if (NS == 0)
			continue;

		for (i = 0; i < N; i++)
		{
			combine_dist (N, L - 1, robustdist + i * L, robustDist, i);
		}

		doprint_stat (N * (L - 1) + 1, robustDist);

		printf ("%3d %4d %lg\n", minlinks, NS2_min, NS);
		for (i = 0; i < L; i++)
		{
			for (j = 0; j < N; j++)
			{
				printf ("%d ", pathway[i][j]);
			}
			printf ("\n");
		}

		/*
		for (i = 0; i <= N * (L - 1); i++)
		{
			if (robustDist[i] > 0.1)
			{
				printf ("c Number of Robust Steps : %3d\tNumber of Solutions : %lg\n", i, robustDist[i]);
			}
		}
		*/
	}
	fclose (fp);
}
